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Abstract 

In the calculation of cross sections for infrared-safe observables in high energy collisions at next- 
to-leading order, one approach is to perform all of the integrations, including the virtual loop 
integration, numerically. One would use a subtraction scheme that removes infrared and collinear 
divergences from the integrand in a style similar to that used for real emission graphs. Then one 
would perform the loop integration by Monte Carlo integration along with the integrations over 
final state momenta. In this paper, we explore how one can perform the numerical integration. We 
study the iV-photon scattering amplitude with a massless electron loop in order to have a case with 
a singular integrand that is not, however, so singular as to require the subtractions. We report 
results for N = 4, N = 5 with left-handed couplings, and N = 6. 
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I. INTRODUCTION 



The calculation of cross sections in the Standard Model and its extensions at next-to- 
leading order (NLO) in perturbation theory inevitably involves computing virtual loop Feyn- 
man diagrams. The standard method for this involves computing the loop integrals analyt- 
ically. Once the one loop amplitude is known analytically, the result can be inserted into a 
calculation of the cross section in which integrals over the momenta of final state particles 
are performed numerically. This is the method that was introduced in Ref. [l[ and is used, 
for example, in the packages MCFM jl] and NL0Jet++ (!]. 

This approach is powerful and has been successfully applied to a number of processes 
of experimental interest. There has been considerable progress [ij in expanding the range 
of processes for which an analytical answer is known. 1 One may hope that the analytical 
approach may develop into a completely automatic way of generating scattering amplitudes 
for a wide class of processes. However, the complexity of the results produced by known 
analytical methods grows rapidly with the number of partons involved in the scattering. For 
this reason, there may be limits to the range of processes for which analytical methods are 
useful. 

One wonders whether a wider range of processes might be amenable to calculation if one 
were, instead, to use numerical integration for the virtual loop integrals. In a calculation of 
a cross section, the numerical integration would be performed along with the integrations 
over the momenta of final state particles, so that there would be a single integration over a 
large number of variables, with the integration performed by Monte Carlo style numerical 
integration. The purely numerical approach will inevitably have its limitations, just as the 
analytical approach does. However the nature of the limitations will be different. For this 
reason, we believe that one should try to develop the numerical method as far as possible 
and see how far back the limitations can be pushed. Eventually, this should involve trying 
several variations on the basic theme of performing the integrations numerically. 

There are already some methods available for doing the virtual loop integrals numerically. 
In one method p, @] , one performs the integral over the energy flowing around the loop an- 
alytically by closing the integration contour in the upper or lower half plane and evaluating 
the residues of the poles in the complex energy plane that come from the propagator de- 
nominators. This is a purely algebraic step. Then the integral over the space momentum 
is performed numerically. There are infrared divergences, but these cancel inside the inte- 
grals between real and virtual graphs that make up a NLO cross section. This method has 
been applied to e + e~ — > 3 jets and is completely practical in that application. For more 
complicated processes, we do not know how to arrange a calculation in this style without 
subtractions. One could add subtractions to the method of Refs. jE], @], but for this paper 
we have chosen a different approach. 

Another method [7| involves transforming the loop integral into the standard Feynman 
parameter representation that one uses for analytically evaluating such integrals. Then the 
integral over the Feynman parameters is to be performed numerically. This method shows 
promise, but is limited by the complexity introduced by expanding the numerator func- 



Here "known" may mean that there exists a computer program to calculate the desired scattering am- 
plitude in terms of known master integrals or other special functions. There are many approaches, some 
of which are conventionally called "semi-numerical" because parts of the calculation involve a numerical 
approach. 
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FIG. 1: Feynman diagram for the iV-photon amplitude. 



tions involved. The method introduced in this paper makes use of the Feynman parameter 
representation while avoiding the complexities introduced by the numerator function. 

This paper represents the second step of a program for calculating virtual loop integrals 
numerically. In the first step [0], we attacked the problem of infrared divergences. Typi- 
cally, the integrals that one wants to evaluate have infrared divergences associated with the 
momenta of particles in the loop becoming collinear with the momenta of massless external 
particles or becoming soft. In Ref. 0, we proposed a subtraction scheme in which one 
subtracts certain counter terms from the integrand, then adds these same counter terms 
back. After summing over graphs and performing the integrals analytically, the counter 
terms that we added back have a simple form that is easily included in the calculation of a 
cross section. Meanwhile, the main integrand minus the counter term integrands combine to 
make an integrand that is free of singularities strong enough to make the integral divergent. 
Thus one can numerically integrate the main integrand minus the counter-term integrands. 

Despite the beauty of this approach, it is one thing to say that one can numerically 
integrate the combined integrand and it is another thing to do it. One needs a practical 
method for doing it. That is what we propose in this paper. 

In order to keep our discussion reasonably simple, we attack a simple problem in which 
the counter terms are not present because the original integral is infrared finite. The problem 
is to compute the amplitude in quantum electrodynamics for scattering of two photons to 
produce N — 2 photons by means of an electron loop. Our formulas include the possibility 
of a non-zero electron mass, but in order to face up to the problem of infrared singularities 
that appear when the electron mass vanishes, we concentrate on the mass zero case. 

The process is illustrated in Fig. [TJ Electron line n in the loop carries momentum I — Q n , 
where Q n is fixed and we integrate over I. The momentum carried out of the graph by 
external photon n is 2 

Pn — Qn+1 — Qn , (1) 



2 Throughout this paper, we adopt a cyclic notation for indices in the range {1, 2, • • • , N}. Thus Eq. {!} 
for n — N is Pn — Qi — Qn- 
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with P% = 0. The propagator denominators provide factors that would lead to logarithmic 
divergences after integration over the soft and collinear regions. However, these divergences 
are cancelled. For each electron line there is a factor (/ — $„)■ Thus the numerator provides 
a factor that removes the soft divergence from the integration region (l — Q n ) 0. Similarly 
at each vertex there is a factor (/ — n+1 ) i n {Pn) (/ — where e n (P n ) is the polarization 
vector of the photon. In the collinear limit (/ — Q n ) — > xP, this gives a factor — x(l — 
x)f n <f. n (P n )f n = — 2x(l — x)f n e n (P n ) ■ P n . This vanishes because e n (P n ) ■ P n = 0. Thus the 
numerator also provides a factor that removes each collinear divergence. The loop integral 
is also finite in the ultraviolet as long as iV > 4. (For N = A the integral is divergent by 
power counting, so a special treatment, discussed later in this paper, is needed.) Thus we 
can present an algorithm that is uncluttered by the counter terms by means of using the 
scattering of two photons to produce N — 2 photons. We reserve the full case of massless 
quantum chromodynamics for a future paper. 



II. THE AMPLITUDE 

We wish to calculate the amplitude for scattering of two photons to produce N — 2 
photons by means of a (massless) electron loop. However, we formulate the problem in a 
more general fashion. The amplitude for any one loop graph can be represented as 

(^c>n (i _ Qi) ,_ m?+i0 - p> 

Here there is a loop with N propagators as illustrated in Fig. [TJ The nth propagator carries 
momentum I — Q n and represents a particle with mass m n . At the nth vertex, momentum 
P n = Qn+i — Qn leaves the graph. In the case to be considered, all the m n and P% vanish, 
but we leave the masses and external leg virtualities open in the general formulas. There 
is a coupling e for each vertex, where e is the charge of the fermion. There is a numerator 
factor that, for the photon scattering case with zero electron mass, has the form 

N(l) = Tr{f N (Ps)(f-Q N )'..{ 1 (P l )(f-Q 1 )} , (3) 

where e,(Pj) is the polarization vector of photon % and e is the electromagnetic coupling. 3 
In other examples, one would have a different numerator function. The only property that 
we really need is that N(l) is a polynomial in I. 

It will prove convenient to modify this by inserting factors im^ in the numerator and the 
denominator, where is an arbitrary parameter that we can take to be of the order of a 
typical dot product Qi ■ Qj. This factor is not absolutely needed for the purposes of this 
paper, but it is quite useful in the case of the subtraction terms to be considered in future 
papers and is at least mildly helpful in the analysis of this paper. With this extra factor, 
we write 

d 4 l iml e N N(l) 1 



J (2tt)4 iml i\ (1 ~ Q l ? ~ ml + iO ' 1 J 



3 Specifically, ei(Pi) for an outgoing photon is e*(Pj, Sj) = e(Pj, — Sj), where Si is the helicity of the photon. 
For an incoming photon, we follow the convention of using a helicity label Si equal to the negative of the 
physical helicity of the photon. Then ej(Pj) is e(— P,, —si). 
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III. REPRESENTATION WITH FEYNMAN PARAMETERS 



In principle, it is possible to perform the integration represented in Eq. fl3j) directly by 
Monte Carlo numerical integration on a suitably deformed integration contour. We have 
looked into this and conclude that it may be a practical method. However, one has to 
pay attention to the singularities on the surfaces (Z — Qi) 2 = m 2 . The geometry of these 
surfaces and of their intersections is somewhat complicated. There is a standard method 
for simplifying the singularity structure: changing to a Feynman parameter representation. 
One way of using this method has been emphasized in the context of numerical integrations 
in Ref. [7|. It is the Feynman parameter method that we explore in this paper. 

The Feynman parameter representation of Eq. (j3J) is 

The denominator here is 

N 

D(l) = £V[(Z - Q t ) 2 - m 2 } + ix°m 2 . (6) 

i=l 

The denominator comes with a "+i0" prescription for avoiding the possibility that the 
integrand has a pole on the integration contour and the +ix Q m 2 ] term serves the same 
purpose. With repeated use of Ylii=i x% = 1 — x°, the denominator can be simplified to 

D(l) = -±-Jl 2 + A 2 (x)} . (7) 



Here 



and 



1 - x° 

N 

1 = ^2^(1 - Qi) (8) 



1 N N 

A 2 ( x ) = -J2 xWSij + i x°x j m 2 , (9) 



2 

t,j=i j=i 



where we have defined 



S l] = {Q l -Q ] ) 2 -m 2 -m 2 . (10) 
We change integration variables from I to I as given in Eq. (jSJ). The inverse relation is 



/ 

With these results, we have 



M = im 2 e N T(N + 1) J^dx 1 ■ ■ ■ j\x N 9 x ' < 1 j (^2 
dH N(l(J,x)) 



x 



(12) 



(2tt) 4 r~ i^ 1 
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Here we understand that the Feynman parameter x° is given by x° = 1 — Yli=i x% ■ 

If we wished to perform the integration analytically, the next step would be carry out 
the integration over I. However for a numerical integration, such a step would be a step in 
the wrong direction. Performing the / integration analytically would require expanding the 
complicated numerator function in powers of I. For this reason, we leave the I integration 
to be carried out numerically after a little simplification. 

The simplification is to change variables from I to a momentum £ that has been scaled 
by a factor A(x) and rotated in the complex plane: 

F(xJ) = |A(x) {(1 - i)^ + (1 + i)P»t) . (13) 

Here k = P^t is the parity transform of L £° = £°, £ j = -£ j for j e {1 ,2,3}. W e have 
defined A(x) for real x and m 2 , — > to be a/A 2 (x) if A 2 (x) is positive and iy— A 2 (a;) if A 2 (x) 
is negative. The square of / is 

I 2 = A\x)tP^ v t . (14) 

Note that £ fl P fll/ £ p is the square of £ with a euclidian inner product and is thus strictly 
positive. 

Our integral now is 



M = -m 2 e N V(N + l) J 



d 4 



N \ /N \ N -3 AJ(1( ... (15) 

N(l(x,£)) 



j^x 1 --- j\x N el^x' <u [f^ 



X 



[A 2 (x)+iO] N - 1 ' 

The function l(x,£) in the numerator function is obtained by combining Eqs. ( flTT) and ( fl3l) : 

1 



F(x,£) 



x° 



N 

|A(x) {(1 - i)F + (1 + i)P^t} + J2 xj Qj 

3=1 



(16) 



It is a somewhat subtle matter to verify that the complex rotations involved in defining £ 
are consistent with the +i0 prescription in the original denominator. We examine this issue 
in Appendix |A] 

Notice that in the numerator function the momentum on line n is 

/"(*, £)-Q*= 3-^0 Il A (*) {C 1 - W + C 1 + + K n( x )] . ( 17 ) 

where 

AT 

We shall meet K n {x) later in Sec. IVII when we study pinch singularities. For the moment, 
we note simply that in the final formula ( fl5l) both the numerator and the denominator are 
invariant under shifts Qi — ► Qi + AQ of the reference momenta Qi. 
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IV. CONTOUR DEFORMATION IN FEYNMAN PARAMETER SPACE 



The integral in Eq. ([15]) is not yet directly suitable for Monte Carlo integration. The 
problem is that the quadratic function A 2 (x) vanishes on a surface in the space of the 
Feynman parameters. Evidently, the integrand is singular on this surface. For x° > 0, A 2 (x) 
does not vanish for real x\ but on the plane x° = 0, A 2 (a;) vanishes for certain real values 
of the other xK If we don't do something about this singularity, the numerical integral will 
diverge. The something that we should do is deform the integration contour in the direction 
indicated by the +i0 prescription. That is, we write the integral as 



N-l 

ZY 



,-^r(Ar + i)/|^ [i+ ^ r+l (19) 

x/V-/V,(Ee<i)det(|) (£• 



Here we integrate over real parameters for i G {0, 1, . . . , A^} with YliLa = 1, so that we 
have displayed the integral as an integral over A" parameters . . . , £ N with £° = 1 — YliLx £*■ 
The integration range is < for i G {0, 1, ... , A^}. The original integral was over real 
parameters x % with Y2iLo xl = 1 with this same range, < x l . The contour is defined by 
specifying complex functions £*(£) for i G {0, 1, . . . , N} with 52 i=0 z % = \. 

In moving the integration contour we make use of the multidimensional version of the 
widely used one dimensional contour integration formula. A simple proof is given in Ref . (fif . 
The essence of the theorem is that we can move the integration contour as long as we start in 
the direction indicated by the +i0 prescription and do not encounter any singularities of the 
integrand along the way. In addition, the boundary surfaces of the contour have to remain 
fixed. Since the surfaces z l = are boundary surfaces of the contour before deformation, 
they should remain boundary surfaces after the deformation. The original integral covers 
the region < x % for i G {1, . . . ,N} and the £ l cover this same range, < £\ Thus we 
demand 

z\0 ^ as f - (20) 

for % G {0,1,..., AT}. 

We adopt a simple ansatz for the contour in the complex z-space: 4 



l + iEf=o^'(0 



m = : \,r s L. • (2D 



Here the rf variables are functions of the integration parameters £\ With this ansatz, the 



4 This is a non-trivial deformation for all £ such that 7y(£) ^ with one exception. If all of the vanish 
except for £™, where then £ ra = 1, then z % = for all % even if rf 1 7^ 0. This possibility does not cause any 
problems. 
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constraint that Y2i ^ = 1 is automatically satisfied: 

N N 

J2e = i =► E / = 1 - ( 22 ) 

i=0 i=0 

In order to satisfy Eq. ( 1201) . we require 

rf(0 ^ as f -> (23) 

for % G {0,1,..., AT}. 

There are certain conditions to be imposed on the contour choice in order to be consistent 
with the "+i0" prescription in the original integral. Note first that 

A' W = AH( V' m) . (24, 

Next, note that A 2 with argument £ + i?y appears in the numerator. In order to analyze 
Eq. (|24p . it is convenient to give a special name S(x) to the quadratic function that forms 
the first part of A 2 in Eq. (jUj), 

TV 

A 2 (:r) =S{x) +i^xV m 2 , (25) 
i=i 

where 

1 " 

5(x) = -J2 ■'■'■' JS o ■ ( 26 ) 

».i=i 

A sufficient condition for the choice of the (£) is as follows. First, we choose 

77° = . (27) 
This is the simplest way to satisfy Eq. ( |23i) for 77°. With this choice for r)°, we have 

N N 

A 2 (£ + i V ) = 5(0 - S( V ) - m ieY,^ + i Y, ^(0 ^(0 + im o - £°) , (28) 

j=l i=l 

where 

«,(0 S ^ = ESW • (29) 
Our condition for the choice of the rf for % e {1, . . . , A^} is that 

AT 

X>^M(£)>o » (3°) 

1=1 

with ^ > except at a point on the boundary of the integration region. 

Suppose, now, that the condition (130]) is satisfied. Do we then have an allowed contour 
deformation? Consider the family of contour deformations A) = Arf (£) with < A < 1. 
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We first consider infinitesimal values of A. We have, to first order in A, 



A 2 (z) 



N 



N 



X 



3=1 



i=l 



/v 



i=i 



+ C»(A 2 



5(0 - Am 2 £° £ ^ + 2A£°(1 - £ Vo E ^' 
i=i 3=1 

N N 

+ iA v\o ™m - 2iA5(o E ^(0 + im o - e°) + <?(a 2 



(31) 



i=i 



In the neighborhood of any point £ with £° > 0, A 2 (z) has a positive imaginary part even 
with A = 0. For £° = 0, the contour deformation gives A 2 (z) a positive imaginary part in 
a neighborhood of any point £ where the real part, 5(0, vanishes. This is the meaning 
of the "+i0" prescription. We may consider that we start with a value of A that is just 
infinitesimally greater than zero, so that the contour does not actually pass through any 
poles of the integrand in the interior of the integration region. 
Now we turn to larger values of A. We have 



A\z) 



A 2 (£ + iA?7(0) 

;i+iA£^(o) 2 



(32) 



Assuming that the ?/(0 are smooth functions, this is a smooth function of £. (Note here 
that 1 + iA^j^O cannot vanish because its real part is 1). Furthermore, when A > 0, 
A 2 (z) is never zero in the interior of the integration region. This is because, according to 
Eq. (I28p . the imaginary part of A 2 (£ + iAr/(0) is positive. Thus l/[A 2 (z)} N ~ l is an analytic 
function of z in the interior of the entire region covered by the family of deformations. For 
the boundary of the integration region, there are some issues of convergence that one should 
check. We do so in Appendix [B] Anticipating the result of this check, we conclude that the 
integral is independent of the amount of deformation and we can set A = 1. 

It is remarkable that the imaginary part of A 2 (z) is not necessarily positive on all of the 
deformed contour. What is crucial is that the deformation starts in the right direction and 
that, as the contour is deformed, it does not cross any poles. 



V. A STANDARD CONTOUR DEFORMATION 

A convenient choice for the deformation function rf (£) for i e {1, . . . , A^} is 

/f(0 = (A/m 2 )f<0 • (33) 

Here A is an adjustable dimensionless constant and m 2 is a parameter with the dimension of 
squared mass that we insert because Sij and thus Wi has dimension of squared mass. Note 
that with this choice the requirement fl23l) that rf(0 vanish when £ l vanishes is automatically 
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met. This deformation gives 

N 

S(Z + VQ) = S(0 - S( V ) + i(A/m 2 ) C MO] 2 • (34) 

i=i 

Evidently the imaginary part of A 2 (£ + iq) has the right sign. 

Eq. (1331) can be thought of as specifying a basic deformation. We can add other defor- 
mations to this. In our numerical work for this paper we have added one more deformation, 
as specified in Appendix 

VI. PINCH SINGULARITIES 

The integrand is singular for £° — > at any real point £ with = 0. We have seen 
in the previous section that the standard contour deformation keeps the contour away from 
this singularity as long as there is some index i G {1, . . . , N} such that > and ^ 0. 
What about a point £ with «S(£) = such that there is no index i G {1, . . . , N} such that 
> and Wi(£) ^ 0. In this case, the integration contour is pinched in the sense that there 
is no allowed contour deformation that can give <S(£ + irf) a positive imaginary part at this 
point £. To see this, recall from Eq. (j28p that, when £° = 0, 

N 

ImS(£ + i77) = £^(0 . (35) 
i=i 

Consider a point £ such that <S(£) = and such that for each index i G {1, . . . , A/"}, 7^ 
implies = 0. For any allowed deformation ry(£), we must have rf = for all i G {1, . . . , A^} 
such that £ l = 0. Thus for each index i G {1, . . . , N}, Wi(£) 7^ implies r] 1 = 0. From 
Eq. (|35p we conclude that Im<S(£ + iry) must vanish at the point in question for any allowed 
choice of the rj l . We conclude that a real point £ with £° = and with <S(£) = is a pinch 
singular point if, and only if, 

CwiiO = for every i G {1, . . . , N} . (36) 

We also note that the point £° = 1, £ l = for i G {1, . . . , A^}, is a pinch singular point. 
This singularity corresponds to the ultraviolet region of the original loop integration. 

With a little algebra, one can translate the condition for a pinch singularity with £° = 0. 
At one of these points, we have 

either f = or Kf - m\ = (37) 

for each i G {1, . . . ,N}, where Kf (£) was given earlier in Eq. (jTSJ). When £° = 0, these 
vectors have the properties that K { — K i+ i = Pi and = 0. Thus Eq. (I3TI) is the well 

known condition for a pinch singularity (see Bjorken and Drell It says that for each 
propagator i around the loop, there is a momentum Ki such that momentum conservation is 
obeyed at the vertices and each Ki around the loop is either on shell or else the corresponding 
is zero and such that the space-time separations Axf = £ l A"f around the loop sum to 
zero. 
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FIG. 2: Illustration of the double parton scattering singularity. 

Notice that the momenta K?(£) appear in the numerator function. According to Eq. (1171) . 
the momentum for line i in the numerator function in the case that £ is at a contour pinch 

(soA(0 = 0)isirf(0. 

There are two types of pinch singular points that are always present if we have massless 
kinematics (with no external momenta collinear to each other) and one more that can be 
present. 

A. Soft singularity 

The first kind of pinch singular point that is always present if we have massless kinematics 
is the one corresponding to a loop propagator momentum that vanishes. If m n = for some 
n then S nn = 0. This means that A 2 (£) = when all of the vanish except for £ n , which 
is then £ n = 1. This is a pinch singular point because all of the z l (^) are fixed: z l = 
for % ^ n, z n = 1. This point corresponds to the momentum of line n in the momentum 
space representation vanishing. In our photon scattering example, there is a singularity at 
this point but because of the zero from the numerator function it is not strong enough to 
produce a divergence. 

B. Collinear singularity 

The second kind of pinch singular point that is always present if we have massless kine- 
matics is the one corresponding to two loop propagator momenta becoming collinear to an 
external momentum. If m n = m n +i = for some n and if P% = (Q n +i — Qn) 2 = 0, then 
S nn = S n +i, n +i = S n +i, n = 0. This means that A 2 (£) = when all of the vanish except 
for £ n and £ n+1 . It also means that w n (£) = u>„+i(0 = 0, so that this is a pinch singular 
point according to the condition (1361) . This point corresponds to the momentum of lines n 
and n + 1 in the momentum space representation being collinear with P n . In our photon 
scattering example, there is a singularity along this line but because of the zero from the 
numerator function it is not strong enough to produce a divergence. 

C. Double parton scattering singularity 

A third type of pinch singular point can be present if a special condition holds for the ex- 
ternal momenta. This singularity corresponds to double parton scattering and is illustrated 
in Fig. [2j Imagine that incoming parton A splits into two collinear partons. Imagine also 
that incoming parton with index B splits into two collinear partons. One of the partons 
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from A and one from B could meet and produce a group of final state partons. The other 
parton from A could meet the other parton from B and produce a second group of final state 
partons. For this to happen, we need at least two external lines in each group of outgoing 
partons produced. Thus we need at least four outgoing external particles. Thus we need 
N >6. 

This picture satisfies the criteria of Eq. (157)) for a pinch singularity. In the Feynman 
parameter space, the singularity occurs along a one dimensional line in the interior of the 
space. We work out where this line is in Appendix [Dl 

Now, the pinch singularity conditions hold only for certain special choices of the external 
momenta. However, if N is large, it is usual that the kinematics is close to a pinch singularity 
condition for some of the graphs. For this reason, in a numerical program, one should check 
for each graph if such a nearly pinched contour occurs. In the event that it does, one should 
put a high density of integration points near the almost singular line. 



VII. ULTRAVIOLET SUBTRACTION 

Some graphs are ultraviolet divergent. For instance, in the photon scattering case, there 
is an ultraviolet divergence for N = 4 (and for N = 2, but we do not consider that case.) 
In the representation (1X91) . the divergence appears as a divergence from the integration over 
Feynman parameters near £° = 1, with all of the other near zero. The reader can check 
that with a numerator function proportional to N powers of the loop momentum, this region 
does give a logarithmic divergence for N = 4. If the graph considered is ultraviolet divergent, 
it needs an ultraviolet subtraction, so that we calculate 



•Mnet = M - M m . (38) 

In a numerical integration, we subtract the integrand of M. m from the integrand of M., then 
integrate. We arrange that the singularities of the integrand cancel to a degree sufficient 
to remove the divergence. The subtraction term is defined in Ref. jsj so that it reproduces 
the result of MS subtraction (which we use because it is gauge invariant). In the photon 
scattering case, the sum over graphs of the subtraction term vanishes, corresponding to the 
fact that there is no elementary four photon vertex. Thus the result after summing over 
graphs does not depend on the MS renormalization scale. 
The subtraction term from Eq. (A. 37) of Ref. (8[ is 

f dH im 2 e 4 [ jVq 1,1,1)- 32 fl^etfQ 32UU 1 ' ^ ) m] 
uv J (2tt) 4 im 2 \ [I 2 - A* 2 e" 4/3 + iO] 4 [I 2 - ^e^l 2 + iO] 4 J 1 j 

Here /i 2 is the MS renormalization scale, which can be anything we like since the net counter- 
term is zero. For the numerator function in the first term, we have adopted the notation 
that the ordinary numerator function N(l) in Eq. (jHJ) is written N(k&, k%, fv2, ki) where k n = 
I — Q n . In this notation, N(l, I, I, I) is the standard numerator function with each propagator 
momentum set equal to /. 

We can now apply the same transformations as for the starting graph to obtain the 
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representation 

d 4 



M uv = -m 2 e 4 T(5) J 



(2tt) 4 [1 + evp^]* 



dx 1 --- J dx 4 9 (^Z %l < *J y 

iV(/, Z, I, I) - 32 n- =1 1 ■ ejiPj) 32 []J =1 T- e^P,) 



Here 



[A 2 4/3 (x) + i0]3 [A 2 3/2 (x)+i0]3 

A\ /3 (x)= -(l- x y^ e -^ + ix (l-x°)ml , 
A 2 3/2 (x) = -(l-x°)Ve- 3/2 + ia; (l-a; )m 2 , 
„o 



(41) 



where we have used x — 1 — Y2j=i x ° ■ m the numerator of the first term, I is a function 

/(•'•• 0- ^' 

/"(*, *) = [|A 4/3 (x) {(1 - i)r + (1 + i)PX}} • (42) 

In the second term, I in the numerator is a function l(x,£), 

H x J) = T ^- [^3dx){(l-i)F + (l + mn] ■ (43) 

The reader can check that for the photon scattering case with AT = 4 the integrand for 
ultraviolet subtraction matches that of the starting graph in the region x° — * 1, so that if 
we subtract the integrand from the counter-term graph from the integrand for the starting 
graph, the resulting integral will be convergent. 



VIII. THE MONTE CARLO INTEGRATION 



We have implemented the integration in Eq. ffl9l) as computer code 10j. The integration 
is performed by the Monte Carlo method. This is a standard method, but it may be good 
to indicate what is involved. First, we note that we do not simply feed the integrand to a 
program that can integrate "any" function. There are many reasons for this, but the most 
important is that we do not have just any function but a function with a known singularity 
structure, a structure that is generic to loop diagrams in quantum field theory with massless 
kinematics. We can take advantage of our knowledge of how the integrand behaves. 

To proceed, we note that we have an integral of the form 

M = fdH /V /V ■ ■ ■ / d ^ 6 \ ~ 1 I f&t) ■ ( 



JO 




In a Monte Carlo integration, we choose N pts points {Ej, £j} at random with a density p(£, £) 
and evaluate the integrand /(£) at these points. Then the integral is 

M= lim J-V 4t44 • (45) 
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The integration error with a finite number of points is proportional to 1/ y^N pts . The coef- 
ficient of 1/ a/ ^points m the error is smallest if 

p(£,0 « const, x • (46) 

That is the ideal, but it is not really possible to achieve this ideal to the degree that one has 
a one part per mill error with one million points. However, one would certainly like to keep 
\f(£, £)\/p(£, £) from being very large. In particular, /(£,£) is singular along certain lines in 
the space of the £ (the collinear singularities) and at certain points (the soft singularities). We 
need to arrange that p is singular at the same places that / is singular, so that f(£, £)/p(£, £) 
is not singular anywhere. Since \f(£, £)| can be very large near other lines associated with 
double parton scattering, we also need to arrange that p(£,£) is similarly large near these 
lines. 

We construct the desired density in the form 



alg 

p(U) = P<WX>jA/(0 • (47) 
j=i 

Here j d 4 £ pe(£) = 1, the sum of the aj is 1, and there are several density functions pj with 
Jd£,Pj{0 = 1- Each pj corresponds to a certain algorithm for choosing a point £. For 
each new integration point, the computer chooses which algorithm to use with probability 
a j. The various sampling algorithms are designed to put points into regions in which the 
denominator is small, based on the coefficients SV,-. We omit describing the details of the 
sampling methods since these are likely to change in future implementations of this style of 
calculation. 

Points I are chosen with a simple distribution pe(£). In the calculation of the numerator, 
we average between the numerator calculated with i and the numerator calculated with —£. 

Having outlined how j\A is calculated by Monte Carlo integration, we pause to suggest 
how the calculation of a cross section (for, say, Higgs production) would work. There one 
would have one-loop amplitudes M. expressed as integrals and one would need to multiply 
M. by a function h(P) of the external momenta that represents a tree amplitude and a 
definition of the observable to be measured. One would need the integral of this over the 
external momenta P. One would perform all of the integrations together. That is, one 
would choose points {P, £, £} and calculate the contributions from the virtual graphs times 
tree graphs to the desired cross section according to 

j= n m J_Y f ty& Pj l h(Pi) • (48) 

Here / is the integrand of M. as above. The function p is the net density of points in £, £, 
and P. Thus what one would use is not M. itself but rather the integrand for M.. 



IX. CHECKS ON THE CALCULATION 



we have implemented the integration in Eq. ffl9|) 
With this code, there are a number of internal checks that can 
be performed on the computation. First, we can replace the real integrand by a function 



As discussed in the previous section, 
as computer code 
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that has soft or collinear singularities but is simple enough to easily integrate analytically. 
Then we can compare the numerical result to the analytical result. This checks that the 
functions Pi(£) and p{£) correspond to the true probabilities with which points £ and £ are 
chosen. Then we can vary the amount of deformation (both for the real integral and for the 
test integrals). When we integrate over a different contour, the integrand is quite different. 
Nevertheless, the (N + 4)-dimensional Cauchy theorem guarantees that the integral should 
be unchanged, provided that the integration is being performed correctly. Thus invariance 
under change of contour is a powerful check. Another check is to change the value of the 
parameter mo- At the start, M. is proportional to m\jm\ and is trivially independent of mo- 
However in the integral as performed, Eq. ffl9|) . mo is deeply embedded into the structure 
of the integrand, so that it is a non-trivial check on the integration that the result does not 
change when we change m . Next, we can replace one of the photon polarization vectors 
e n (P n ) by P n . This gives a non-zero result for each Feynman graph, but should give zero for 
the complete amplitude summed over graphs. Additionally, we can change the the definition 
of the polarization vectors e n (P n ). For reasons of good numerical convergence, we normally 
use polarization vectors appropriate for Coulomb gauge, but we can switch to a null-plane 
gauge. The two amplitudes should differ by a phase, so that \M\ is unchanged. Another 
check is obtained by replacing the vector current by a left-handed current or a right-handed 
current. For even N, the left-handed and right-handed results should be the same, while for 
odd N they should be opposite. For another check, we can reformulate the integral so that 
we do not define i with a scale A(x) _1 . In this formulation, the denominator is 

[S{x) + ix°(l - x°)mo(l + ^P^t)] N+1 . (49) 

The structure of the integral is quite different, but the result should be the same. For four 
photons, there is one additional test: the result should be independent of the renormalization 
parameter //. 

We have subjected the code to these checks. The numerical precision of the result is 
often not high and we have not used every check for every choice of N and external momenta 
and polarizations. Nevertheless, we have found that where we have tried them, the various 
checks are always passed. We note that a better check would be to obtain the same results 
with completely independently written code. We have not done that. 

X. RESULTS 

In this section, we use this code to test how well the method described works. The 
result for a given choice of helicities of the photons has a phase that depends on the precise 
definition of the photon polarization vectors e«. However, the absolute value of the scattering 
amplitude M. is independent of this conventions used to define the e^, so we concentrate on 
\M\. Our convention for defining M. is specified in Eq. (j2J). Since \M\ is proportional to 
ol n I 2 and has mass dimension 4 — N, we exhibit \Ai\ x (i/s^^/a^ 2 in our plots. We 
specify helicities in the form hi, hi, ^3, • • • , h^, where 1 and 2 are the incoming particles 
and, following convention, hi and hi are actually the negative of the physical helicities of 
the incoming photons. 
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FIG. 3: Four photon amplitude. We plot |.M|/a 2 for helicities -H versus 6 — ir/2, where 6 is 

the scattering angle. The curve is the analytical result from Ref. [ill ]. The points are the result of 
numerical integration. 



We begin with N = 4, light-by-light scattering. Here we use the subtraction for the 
ultraviolet divergence in each graph as described in Sec. IVHI For the N = 4 case the result 
is known and has been presented in a convenient form in Ref. [llj. For the two helicity 
choices +++-|- and +-H — , |A^|/a 2 = 8. Our numerical results agree with this. For the 

choice -H , the result depends on the value of the scattering angle 9. In Fig. [3]we exhibit 

the prediction of Ref. [ll[ versus 9 — it/ 2 as a curve and a selection of points obtained by 
numerical integration as points with error bars. The error bars represent the statistical 
uncertainty in the Monte Carlo integration. It is not easy to see the error bars in the figure. 
The fractional errors range from 0.0022 to 0.0034. The points were generated using 10 6 
Monte Carlo points for each of six graphs. 



We turn next to N = 5. Since the five photon matrix element vanishes after summing 
over graphs, we use a massless vector boson that couples to the electron with the left-handed 
part of the photon coupling. The final state phase space has four dimensions, which does 
not lend itself to a simple plot. Accordingly we have chosen an arbitrary point for the final 



A. 



N = 4 



B. N = 5 
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FIG. 4: Five vector boson amplitude. We plot y/s \M\/a 5 ^ 2 with helicities H h+. The vector 

boson is massless and couples to the electron with the left-handed part of the photon coupling. An 
arbitrarily chosen final state was rotated about the y-axis through angle 9. 

state momenta {^3,^4,^5}: 

p 3 = (33.5,15.9,25.0) , 

p 4 = (-12.5,15.3,0.3) , (50) 
p 5 = (-21.0,-31.2,-25.3) . 

We take photon 1 to have momentum pi along the — z-axis (so the physical incoming mo- 
mentum is along the +z-axis), and we take f>2 along the +z-axis. Then we create new 
momentum configurations by rotating the final state through angle 9 about the y-axis. In 
Fig. HI we plot computed values of y/s \ A4\/a 5 ^ 2 versus 9. The points were generated using 
10 6 Monte Carlo points for each of 24 graphs. 

C. N = 6 

Finally, we compute the six photon amplitude. Here analytic results are known for the 
helicity choices ++++++ and ++++H — : for these helicity choices, the amplitude should 

vanish [l2j]. There is also a non-zero analytical result for the choice +H [l3[ • We 

compute s \A4\/a 3 for these helicity choices and also for H hH — , for which we know of 

no analytic result. Following what we did for iV = 5, we choose an arbitrary point for the 
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FIG. 5: Six photon amplitude. We plot s|A4|/a 3 . An arbitrarily chosen final state was rotated 
about the y-axis through angle 9. The points are the result of numerical integration. At the top, 

the points in the range 10000 to 25000 are for helicities +H and are compared with the 

analytical results of Ref. 13j |. In the middle, the points in the range from 2000 to 8000 are for 

helicities H H — • There is no analytical result for this helicity combination. At the bottom, 

we show numerical results for helicities ++++++ and ++++H — . According to Ref. [121], the 
amplitude should vanish for these helicity choices. The results for ++++++ are computed at 
9 = 0, 0.2, 0.4, . . . while the results for ++++H — are computed at 9 = 0.1, 0.3, 0.5, 



final state momenta {p3,P4,P5,Pe}'- 

p 3 = (33.5,15.9,25.0) , 
p 4 = (-12.5, 15.3, 0.3) , 

p 5 = (-10.0,-18.0,-3.3) , [ ' 

p 6 = (-11.0,-13.2,-22.0) . 

We choose p\ and f>2 as we did for N = 5. Then we create new momentum configurations 
by rotating the final state through angle 6 about the y-axis. In Fig. [3, we plot computed 
values of s \M.\/cc' versus 6. For ++++++ helicities, we compute the amplitude at 6 = 
0, 0.2, 0.4, .... The results are consistent with the known result of zero. For ++++H — 
helicities, we compute the amplitude at 9 = 0.1, 0.3, 0.5, .... The results are again consistent 

with zero. For +H we compare the numerical results to the analytical results of 

Ref. [13j (top curve) and find good agreement. For the helicity choice H hH — , the results 

lie in the range from 2000 to 8000 and exhibit some variation as the final state momenta 
are varied. We do not have an analytical curve with which to compare. The points were 
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generated using 10 6 Monte Carlo points for each of 120 graphs. 5 



XI. CONCLUSIONS 

In the calculation of cross sections for infrared- safe observables in high energy collisions at 
next-to-leading order, one must treat the real emission of one parton beyond the Born level 
and one must include virtual loop corrections to the Born graphs. Most calculations follow 
the method of Ref. [l[, in which the integration over real emission momenta is performed 
numerically while the integration over virtual loop momenta is performed analytically. One 
can, however, perform all of the integrations numerically. 6 

In one approach to the calculation of loop diagrams by numerical integration, one would 
use a subtraction scheme [s| that removes infrared and collinear divergences from the in- 
tegrand in a style similar to that used for real emission graphs. Then one would perform 
the loop integration by Monte Carlo integration along with the integrations over final state 
momenta. In this paper, we have explored how one would perform the numerical integration. 
We have studied the iV-photon scattering amplitude with a massless electron loop in order 
to have a case with a singular integrand that is not, however, so singular as to require the 
subtractions of Ref. [sj . 

One could perform the integration either directly as an integral J d 4 l or with the help a 
different representation of the integral. We have chosen to explore the use of the Feynman 
parameter representation because it makes the denominator simple. We have found that this 
method works for the cases of 4, 5, or 6 external legs. There is, in principle, no limitation to 
the number of external legs. However, for more external legs, the integrand becomes more 
singular because the denominator is raised to a high power, N — 1. This is evident in our 
results by examining the growth of the integration error as N increases. 

In many practical calculations, the partons in the loop can have non-zero masses and 
the partons entering the loop can be off-shell. These possibilities make the analytical re- 
sults more complicated, but we expect that they make the numerical result more stable by 
softening the singularities. 7 However, we leave exploration of this issue for later work. 

It is remarkable that the method presented here works for quite a large number of external 
legs. However, we expect that the method can be improved. One approach lies in making 
a sequence of small improvements that together amount to a big improvement. Along 
these lines, one can work on the algorithm for deforming the integration contour and on the 
sampling methods used for choosing integration points (which methods we have not discussed 
here). Alternatively, one can look for a different representation of M. as an integral. One 
could use an integral transformation other than that provided by the Feynman parameter 
representation or one could use a more direct representation of the integral. In particular, 
the representation of Refs. 0, @] recommends itself. Here we turn the integral J d A l into a 
three dimensional integral J dH <5+(/ 2 ) that is rather similar to what one has for real parton 



5 This takes a bit under an hour for each point on one chip of our computer, but we note that computer 
timings are dependent on the computer and the compiler. 

6 That it is practical to do so has been demonstrated for the case of three jet production in electron-positron 
annihilation p, 0|. However, the method used there does not extend well to the hadron-hadron case. 

7 Having an unstable massive particle as an incoming parton would be an exception, since this can put new 
kinds of singularities into the integrand. 
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emissions. In contrast to Refs. [5|, |6j, however, one would use explicit subtractions. 8 We 
expect that this method or something similar might be superior to the Feynman parameter 
method used in this paper because for large ./V one does not raise a denominator to a high 
power. 

The challenge would be to find a representation that is simple and for which the integrand 
is well enough behaved that one can get numerical results for, say, N = 12. We hope that 
others might accept this challenge. 
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APPENDIX A: WICK ROTATION 

In this appendix, we explain the contour deformation necessary to obtain the scaled and 
rotated momenta of Eq. (fT3]) . The simplest procedure is to start by rotating the space-parts 
of the vector I, 

1° = k° , V = e- w k J k = 1, 2, 3 . (Al) 

Here the components of k are real. We start with 9 = and increase 9 until 9 = ir/2. Thus 
we rotate the / contour. Throughout the rotation, I 2 has a positive imaginary part. At the 
end, I 2 + A 2 (x) becomes k' 1 P flu k l/ + A 2 (x), where k^P^k" is the euclidean square of k. 

The next step is to rotate all of the components of k by half of the phase of A 2 (x), so 
that after the rotation, k^P lJLV k u has the same phase as A 2 (x) (which itself has a positive 
imaginary part). 

Finally we rescale the components of k by the absolute value of A 2 (x). Thus 

k" = A{x)F . (A2) 

The net transformation is that of Eq. (fT3|) . At all stages, the imaginary part of the denom- 
inator is positive. 



APPENDIX B: CONTOUR DEFORMATION 



In this appendix, we exhibit some details of the argument that the integration over 
Feynman parameters is left invariant by the contour deformation. We start with 



Jd^detAF(z^X)) . 



(Bl) 



We understand that S. Catani, T. Gleisberg, F. Krauss, G. Rodrigo, and J. Winter are working along 
these lines 14 1. 
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Here z(£, A) specifies the deformed contour, 



/pi pi JV 

Jo Jo =1 



(B2) 



the matrix A is 

A ) = • < B3 > 

and is the integrand, 



In order to take care with what happens at the integration boundary, we define a function 
R(£) that measures how far the point £ is from the boundary, 

/?(0 = min(V,...,£V-|]ej • (B5) 

The boundary is at R(£) = 0, and in general < R(£) < l/(N + 1). Then 

I = lim/(r) , (B6) 

where 

J(r) = Jdi 9(R(0 > r) det A F(z(£, A)) . (B7) 

Now let us make a small change 5X in the deformation parameter. If we can prove that 
the corresponding change 51 of the integral vanishes, then the integral is the same for A = 1 
as it was for an infinitesimal A. We calculate 5I(r) for non-zero r. As shown in Ref. 
5I(r) is the integral of a total derivative. Thus we get an integral over the boundary of the 
contour, 

SI(r) = Jdi det A S(R(C) - r) 6R(£) F{z($) . (B8) 

Here 

where B is the inverse matrix to A, 

E4^ = 4- (bio) 

3 

We think of A as producing a vector 5z from a vector 5£, 8z l = J?. A* 5£ J : . Then we can 
think of i? as producing a vector 5£ from the vector 5z given by the change of z under the 
change of deformation, 
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This justifies the name SR for the combination 



6r = T,t£ 6 ? ■ ^ 



Given the ansatz (j2T|) for z(£), the variation 5z takes the form 

i(5A 



[l + iAEf=i^(01 5 



A' 



(B13) 



We build into the definition of the contour deformation the requirement that as any £ fe 
vanishes, the corresponding r] k also vanishes, with r\ k oc £ fc . Then Sz k oc £ k in this limit. 
Also, when 1 - ^£ fc -> 0, it follows from Eq. flBT3l) that ^ 5z fc oc 1 - The result 

is that as we approach a boundary of the integration region, R(£) — > 0, the function 5R(£) 
vanishes, with 

5R(0 = R(0 x , (B14) 

where h(£) is non-singular. Thus 

81{r) = rx J d£ det A 5(R(£) - r) h(£) F{z{£)) . (B15) 

The factor r would seem to imply that 5I(r) — > as r — > 0. However, we should be careful 
because F(z(£)) is singular near the boundary of the integration region. To examine this 
issue, we note that 

,1/(7V+1) 

I = dr J(r) , (B16) 



o 



where 

7(r)= / d£ det A 5(R(£) — r) . (B17) 



Were it not for the numerator function (and the UV subtraction in the case N = 4), the 
integral for / would be logarithmically divergent. Generically, a one loop integral could 
produce two logs, so that I(r) would have a singularity log(r)/r for r — > 0. However, 
the numerator factor (and UV subtraction if needed) produces an extra factor of r. Thus 
I(r) oc r° log (r) for r — > for some if. The power counting for 5I(r)/r, with its extra non- 
singular factor /t(£), is the same. Thus 5I(r) is proportional to r times possible logarithms 
of r as r — > 0. 

We conclude that when we make an infinitesimal change of contour with the properties 
specified in this paper, the variation of the integral vanishes, 

51= lim 5i(r) = . (B18) 

i — >o 

Thus the integral on the deformed contour is the same as on the original infinitesimally 
deformed contour. 
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APPENDIX C: EXTRA DEFORMATION 



Here we return to the question of contour deformation. We study a problem that can 
occur with the standard deformation. We start by stating the problem rather abstractly. 
Let C be a subset of {1,2,..., N} and let B be its complement. Suppose that 

Sij = for i,jeB . (CI) 

We consider the following limit. Define 



and let 



? = tc& for j EC 



(C2) 



(C3) 

e = (i-e-&)e B for jeB , 

so that 

E^E^ 1 • (° 4 ) 

jec jeB 

Then we consider the limit £c —* with £° = 0. Thus the £ l for i £ B are big and the £ J for 
i E C are little. 

In the limit <^£ — > 0, <S(£) becomes 

5 (o = E E ^ + 1 E = & E + > (G5) 

where 

• (C6) 

jeB 

If we adopt the standard contour definition from Eq. (1331) . we have 



S(£ + in) = Zc E + i A & E &[^(&)1 2 + • ( C7 ) 

We see that the surface £c = with £° = is a singular surface of the integrand. In fact, 
it is a pinch singular surface. For a generic point £, <S(£ + ir/) vanishes linearly with ^£ as 
6c -0. 

This generic behavior is fine from a numerical point of view. However, we would like to 
avoid having <S(£ + irj) vanish faster than linearly as £c 0. The real part of S(£ + irj) can 
easily vanish quadratically as <^£ — > 0. The components Sy can have either sign, so that for 
some points £ the particular linear combination Y2iec ££^*(6b) can vanish. It is harder for 
the linear contribution to the imaginary part of <S(£ + it]) to vanish. However, if the set B 
has more than one element, it is possible for u>i(£e) to vanish for some particular index % — I 
at some particular value of £g. Then if all of the £ £ vanish except for i = I, we will have 



E&[^)] 2 = Sm^)] 2 = • ( Q 
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For this choice of the £ l , we will have A 2 (£ +ir/) = 0(£c) if we take the standard deformation. 
One might think that having an esoteric integration region in which the integrand is extra 
singular is not a problem. However, in a numerical integration it is a problem. One possibility 
is to put extra integration points in the region of extra singularity, but a more attractive 
possibility is to fix the contour deformation so as to better keep the integration contour away 
from the singularity. At the same time, we need to avoid letting the jacobian det(dz/d£) in 
Eq. ffl9|) become singular. This is the strategy we will pursue. 

Until now, we have followed a rather abstract formulation of the problem for the reason 
that the same abstract problem occurs in several ways in the subtraction terms defined in 
Ref. jl] to take care of infrared divergent graphs. In this paper, however, we are concerned 
with infrared finite graphs representing photon scattering with a massless electron loop. The 
problem is associated with the region in the original loop integral in which lines n and n + 1 
are nearly collinear. With massless kinematics, S nn = SVt+i,n+i = Sn,n+i = 0- Thus the 
matrix has the special form with B = {n,n + 1} and £ consisting of all index values 
% G {1, . . . , N} other than n and n + 1. 

We will seek a supplementary deformation r? n (£) and ?? n+1 (£) that we can add to the 
standard deformation. In the following, we consider n to be any fixed index value in the 
range 1 < n < N. There will be an analogous deformation for any n. For reasons that will 
become apparent, we want to have just one of these added deformations fj(£) for any value 
of £. For this reason, we will arrange that r/ n (£) and f/ n+1 (£) are nonzero only in the region 

n n -. c + c +1 > h € > C + \ C +1 > C' 1 ■ (C9) 

It is easy to verify that the various regions TZ n are non-overlapping. 

Given that there is an added deformation and f/ n+1 (£), there is an added contribu- 

tion to Im^ that has the form 



AIm<S(£ + i V ) = {V ns n,i + V n+1 Sn+i,i} 



(CIO) 

iec 



When we include the standard contour definition from Eq. fl33|) . the total Im<S is 

lmS(£ + irj) =Y,C {v n S n ,i + V n+1 S n+ i, + (A/m 2 )M£)] 2 } 

iec , cn . 



ieB 

,2\ 



We take the extra deformations to be of the form ±(A/m ) 

rf = (A/m 2 ) g(0 , 
V n+1 = -(X/m 2 )g(0. 



„,.._2,.,~ ( C12 ) 



Here g(£) is a function to be defined. 

With the definition ( 1C12I) for the added deformation together with the standard defor- 



mation, we have 

im<s(£ + ir,) = (\/m 2 )J2c {{Sn,i - s n+ iM0 + k(0] 2 } 

_ ieC (C13) 
+ ^(A/m 2 )f[^(0] 2 • 



iec 
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The second term here is always positive but vanishes quadratically with £c in the limit 
£c — > 0, so it is too small to help us in this limit. In the first term, the parts proportional 
to [u>j(£)] 2 are always positive and for typical values of £ n and £ n+1 vanish linearly with £c. 
However, u>j(£) for some value of i can vanish in this limit for a particular value of £ n and 
^n+i This is the reason for adding the new deformation specified by 
We need to ensure that 

(s n>l - s n+h M0 + MOF > o (ci4) 

for all i G C and for all £. Here we really want a ">" relation and not just a ">" relation, 
which could be satisfied with = 0. 

This goal can be accomplished by a straightforward construction. First, we need some 
notation indicating certain sets of indices. Let C + be set of indices in C such that S n +i,i > S n> i 
and let £_ be set of indices in £ such that S n> i > SVh-i,*- The union of C + and £_ is all of 
C (Here we assume that the external momenta do not lie on the surface S n+ i^ = S n ,i for 
any % in C) 



Define 



9+iO = m i n 



9-(0 = mm- . 

Then Eq. ( 1C14I) requires that 

- g40 < 9(0 < ■ (C16) 

Of particular interest is the requirement for a special point for which = 0. If i G 

then g~(£) = at this point and the requirement is that g(£) be positive (but not too 
positive). If i G £+, then g+(£) = at this point and the requirement is that g(£) be 
negative (but not too negative). 

These restrictions are not very restrictive in the case that either g+(£) or <?-(£) is large or 
in the case that one of the sets C± is empty (in which case we interpret the corresponding 
g± to be infinite). In order to ensure that the deformation not be too large, we can also 
impose 

- X g < (X/m 2 ) g(0 < X 9 ■ (C17) 
where X g is a parameter that could be chosen to be A. Defining 

~g±(0=mm[g±(0,m 2 X g /X} , (C18) 

our requirement is 

-g-(0<9(0<9+(£) ■ (C19) 
It is easy to satisfy Eq. (10190 . We set 

g(0 = H(0C 9 [g + (0-g-(0] , (C20) 
where C g is a parameter in the range < C g < 1 (possibly 1/2) and 

H($ = (2c + 2c +1 - 1) 4(c - r +2 xc +1 - 

x 6 (e + C +1 > \) o(c > C +2 ) 0(C +1 > C 



(C21) 
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The purpose of H(£) is to restrict the range of £ for which g>(£) 7^ to the desired region 
lZ n , Eq. (1C9I) . There is also a factor that becomes 4£ n £ n+1 in the limit £c — > 0. This factor 
turns off the deformation as £ n — > or £ n+1 — > 0. Notice that 

< < 1 . (C22) 

With the use of this property, it is evident that the definition (1C20j) satisfies Eq. (1C19I) . 
Suppose that that there is an index % G C + and an index j G £_ such that 

*Sy r)-t-1 ^ Sj 77 ^ , , „ 

S >0 5 <0 (C23) 

and that 

<s 2 . (C24) 

Then there is approximately an effective contour pinch, in the sense that both g + (£) and 
<7~(0 are close to vanishing when all of the £ fc for k G £ are very small and 



q _ q (C25) 

tn+l _ ~7> gyt 

Q _Q_|_Q_Q 
*- , i,n+l u i,n > "j,*! ,J j,n+l 

The contour is pinched already along the whole collinear singularity line £ fc = for k G C, 
but this is an extra pinch that prevents the deformation of this section from being effective. 
For this reason, one should put extra integration points in the region near this point. 

It is instructive to examine the functions u>j(£) for i G £ in the limit that £° and all of 
the £ J for j 6 £ vanish (and assuming massless kinematics). Then the only two x l that 
are non-zero are £ n and £ n+1 = 1 — £ n . Following the notation of Eq. (1C6j) . we can call the 
limiting function Wi(^ n ). We would like to know for what value of £ n (if any) this function 
vanishes. We have 

MC) = s in e + SW(1 - C) • (C26) 

Evidently w«(£ n ) will vanish for some £ n in the range < £ n < 1 if and only if Si n and S^n+i 
are non-zero and have opposite signs. 

Thus we need to know something about the signs of S in and S ijn+ i. First, we note that 
if % — n — 1 then S{ n = 0. Furthermore, if all of the particles i, % + 1, . . . , n — 1 are final state 
particles, then 

fn-\ \ 2 

(C27) 

If two of the particles i, i + 1, . . . , n — 1 are the two initial state particles, then 

Sin= fi>^ >0 • (C28) 




If exactly one of the particles i, % + 1, . . . , n — 1 is an initial state particle, then 5j n < 0. The 
proof of this amounts to showing that if a massless particle A turns into a massive particle 
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A' by exchanging a momentum Q and a massless particle B turns into a massive particle B' 
by absorbing the momentum Q, then Q 2 < 0. We omit the details. 

Given these results and the analogous results for Si jTl+ i we can conclude that Si n and 
Si tU+ i are non-zero and have opposite signs when i is neither of n — 1 or n + 2 and external 
photon n is an incoming particle. 

Supposing that photon n is an incoming particle, then the propagator index i is in 
with Si n > and SVn+i < if all of the external particles i, i + 1, . . . , n — 1 are final state 
particles. If, on the other hand, one of them is the other initial state particle, then Si n < 
and S^n+i > and i is in C + . 

As we have seen, when photon n is an incoming particle, the zero of u»i(£ n ), 



££) = c Sl ' n+ \ > ( C 29) 



lies in the integration range, < £ n < 1. In this case, we can say more about the location 
of this zero. Let the index of the other incoming photon be n' . Define 

n'-l 

r = -2P n , ■ Pk/s . (C30) 

j=n+l 

Then if the index i is in the range n' + 1, . . . , n — 1, we have S^ n > and S itn+ i < so 
i 6 C-. Then one can show that 

> r . (C31) 

On the other hand, if the index i is in the range n+2, . . . , n r , we have Si >n < and Si >n +i > 
so i 6 C + . Then one can show that 

< r . (C32) 
To prove Eq. (1C32I) . write each outgoing momentum in the form 

P { = - ai P n - biP n > + Pj (C33) 

where Pj ■ P n = Pj ■ P n > = and < a, < 1 and < \ < 1. Then for i in the range 
n + 2, . . . , n' we have 

= ( E ^) = f E «*) f E s + f E f) ■ < C34 ) 

\i=n+l / \j=n+l J \j=n+l / \j=n+l / 

For Si n we & dd one more particle n with a n — —1, b n — 0, and no Pj. Thus 

*n = (^p] = ~ fl - £ %) ( £ 6,) S +(fpj) . (C35) 

\j'=n / V i=n+l / \j=n+l J \j=n+l J 

Then 

\ 2 



^ = ^ — — • (( ' ;J)(,) 
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The {P T ) 2 term in the numerator is negative. Thus 

n— 1 n— 1 

£w<Z>i<E H ;- T • (C37) 

The proof of Eq. (IC31I) is similar. 

Thus in the limit that all of the £ 3 for j e £ are very small, the qualitative nature of the 
deformation constructed here is quite simple. We deform £ n into the upper half plane for 
£ > r and into the lower half plane for C, < r. 



APPENDIX D: DOUBLE PARTON SCATTERING SINGULARITY 

As discussed briefly in Sec. IVI C\ a pinch singular point corresponding to double parton 
scattering can be present if a special condition holds for the external momenta. This singu- 
larity is illustrated in Fig. [2j Imagine that incoming parton with index A carries momentum 
—Pa such that P\ — and that parton A splits into two collinear partons with labels A 
and A + 1. That is —Ka{Q = —(1 — xa)Pa and Ka+i(£) = —xaPa- Imagine also that 
incoming parton with index B carries momentum —Pb such that Pj = and that parton 
B splits into two collinear partons with labels B and B + 1. That is —Kb{£) = (1 — xb)Pb 
and Kb+i(£) = —xbPb- (Here the x's are momentum fractions, not Feynman parameters.) 
Partons A + 1 and B could meet and produce a group of final state partons with labels i 
in a set A = {A + 1, . . . , B — 1}. Partons B + 1 and A could meet and produce a group of 
final state partons with labels i in a set B = {B + 1, . . . , A — 1}. Thus 



^P i= - x A P A - (1 - x B )P B , 
- (1 - x a )Pa ~ x b Pb ■ 



It is convenient to write Eq. (IDIi in terms of the internal line momenta Qi using Eq. (pQ). 
We note immediately that ^2 iej ^Pi = Qb — Qa+i and ^ igB P = — Qb+i are timelike 
vectors. Thus 

Q + ' ^ n (D2) 

Notice that for this kind of singularity to occur, we need at least two external lines in set 
A and two in set B. Thus we need at least four outgoing external particles. Thus we need 
N > 6. 

Given the external momenta, the momentum fractions xa and xb are determined. When 
rewritten in terms of the Q i: the first of Eq. (IDip reads 

Qa+i-Qb+i=xa(Qa+i-Qa)-xb(Qb+i-Qb) , (D3) 



while the second equation in fIDlj) is just the negative of this. Take the inner product of this 
with (Qb+i — Qb) and use 

2(Qa+i - Qa) ■ (Qb+i -Qb) = S , (D4) 
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where 

S = Sa,b+i + Sa+i,b — Sa,b — £>a+i,b+i ■ (D5) 
Also note that (Qb+i — Qb) 2 = and 

2(Qa+i — Qb+i) • {Qb+i — Qb) = Sa+i,b — Sa+i,b+i ■ (D6) 
These relations give 

Sa+i,b — Sa+i,b+i m „ N 

x A = = • (D7) 



We similarly derive 



Sa,b+i — Sa+i,b+i mc \ 
x B = = • (D8) 



The kinematic conditions require that 

£^ = 0, (D9) 

i&A 

where P? is the part of Pi transverse to Pa and Pb- (In this frame the sum of the transverse 
momenta of all of the final state particles vanishes, so the sum of the P^ T for the particles 
in set B also vanishes if the sum for set A vanishes.) The condition for this to happen is 
obtained by squaring both sides of Eq. (ID3I) and inserting the solutions for xa and xb- This 
gives 

Sa+i,b Sa,b+i — Sa,b Sa+i,b+i = . (D10) 

That is, 

det( Q SA ' B q Sa+1 ' b )=0 . (Dll) 

\^A,B+l &A+l,B+l J 

Recall that in order to have a double parton scattering singularity, Sa+i,b > and 
Sa,b+i > 0. The determinant condition then implies that Sa,b an d Sa+i,b+i have the 
same sign. In fact, this sign must be negative. To see this, one may note that, because of 
Eq. (IDlOj) . two alternative expressions for xa are also valid: 



_ S a +i,b _ —Sa+i,b+i _ Sa+i,b — Sa+i,b+i m-in\ 
Xa - 7; ^ - 7; ^ - = • 

&A+1,B — dA,B £>A,B+1 — OA+l,B+l o 

Using the first of these, we see that xa > implies that Sa+i,b ~ Sa,b > 0. But then 
xa < 1 implies that Sa,b < 0. We conclude that for a double parton scattering singularity, 
Sa+i,b > and S a +i,b > 0, S a ,b < and S a +i,b+\ < 0. 

What does this mean in terms of solving Eq. fl36l) ? We demand that Eq. fl36l) hold for 
nonzero £ , (, A+1 , (, B and £ B+1 with all of the other £ l = 0. Thus we need Wb = wb+i — 0, 
or 

Sa,b Sa+i,b ) ( £ A 
Sa,b+i Sa+i,b+x) \i A+1 



. (D13) 



Similarly we need wa = wa+i = 0, or 



Sa,b Sa,b+i \ ( £ B 

£>A+1,B Sa+^B+I J \£, B+1 



(D14) 



29 



One can solve these if the determinant of the matrix is zero, that is if Eq. (1D10I) holds. If it 
does, the solution with £ A + £ A+1 + £ B + £ B+1 = 1 is 



where 



i A = fAX , 

i A+l = {l-f A )x , 

e = f B {i-x) , 

^ +1 = (1-/B)(1-X) 



n Sa+i.b — Sa+i,b+i Sa+i,b — £>A+l,B+l 

J A - 



(D15) 



B 



Sa+i,b ~ Sa,b Sa,b+i — Sa+i,b+i S (r\-\a\ 

£>A,B+1 _ —JA+1,B+1 _ >JA,B+1 — JA+1,B+1 

Sa ,b+i — Sa,b Sa+i,b — Sa+1,B+1 S 



That is, J'a = xa and fg — Xb- In order for the pinch singularity to be inside the integration 
region, £ , (, A+1 , (, B and £ B+1 need to be positive. Thus we need to choose x in the range 

< x < 1 . (D17) 

It is of interest to work out the momenta Eq. (fl8l) . when the external momenta obey 

the condition (1D3I) for a double parton scattering singularity and the Feynman parameters 
£ are given by Eq. (ID15I) . One finds 



K A (0 = (i - f A )P A 
K A+ i(0= -IaPa , 
K B (0 = (l-f B )P B 

K B+ i(0= -SbPb ■ 



(D18) 



These are, of course, the relations we started with. 

We learn that if the determinant condition (1D10I) and certain sign conditions hold, A 2 (£) 
has a pinch singularity along a line that runs through the middle of the integration region. 
Now, the pinch singularity conditions hold only for certain special choices of the external 
momenta. However, one can easily be near to having a pinch singularity. For this reason, in a 

2 

numerical program, one should check for each graph if \Sa+i,b Sa,b+i~ Sa,b Sa+i,b+i \ S , 
with the required sign conditions, for some choice of indices A and B. In that event, one 
should put a high density of integration points near the "almost" singular line. 
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